Thermal decomposition of a honeycomb-network sheet — A Molecular Dynamics 

simulation study 
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The thermal degradation of a graphene-hke two-dimensional honeycomb membrane with bonds 
undergoing temperature-induced scission is studied by means of Molecular Dynamics simulation 
using Langevin thermostat. We demonstrate that at lower temperature the probability distribution 
of breaking bonds is highly peaked at the rim of the membrane sheet whereas at higher temperature 
bonds break at random everywhere in the hexagonal flake. The mean breakage time r is found 
to decrease with the total number of network nodes A*' by a power law r oc N~°'^ and reveals an 
Arrhenian dependence on temperature T. Scission times are themselves exponentially distributed. 
The fragmentation kinetics of the average number of clusters can be described by first-order chemical 
reactions between network nodes Ui of different coordination. The distribution of fragments sizes 
evolves with time elapsed from initially a (5-function through a bimodal one into a single-peaked again 
at late times. Our simulation results are complemented by a set of l''*-order kinetic differential 
equations for rii which can be solved exactly and compared to data derived from the computer 
experiment, providing deeper insight into the thermolysis mechanism. 
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I. INTRODUCTION 

Thermal degradation and stabilization of polymer sys- 
tems has been a long-standing focus of research from 
both practical and fundamental viewpoints Plastic 
waste disposal has grown rapidly to ecological menace 
prompting researchers to investigate plastic recycling by 
degradation as an alternative On the other hand, 
degradation of polymers and other high molecular weight 
materials in different environments is usually a major 
limiting factor in their application. Thermal degrada- 
tion (or, thermolysis) plays a decisive role in the design 
of flame-resistant polyethylene and other plastic materi- 
als [1]. Another interesting aspects for applications in- 
clude reversible polymer networks [1, and most no- 
tably, graphene, as a "material of the future" that shows 
unusual thermomechanical properties 0, Recently, 
with the rapidly growing perspective of exploiting bio- 
polymers as functional materials 0, |^, the stability of 
such materials has become an issue of primary concern 
[13, [m as, e.g., that of double-stranded polymer decom- 
position 

Most theoretical investigations of polymer degradation 
have focused on determining the rate of change of aver- 
age molecular weight |13l423 |. The main assumptions of 
the theory are that each link in a long chain molecule 
has equal strength and equal accessibility, that they are 
broken at random, and that the probability of rupture is 
proportional to the number of links present. Therefore, 
all of the afore-mentioned studies investigate exclusively 
the way in which the distribution of bond rupture prob- 
ability along the polymer backbone affects the fragmen- 
tation kinetics and the distribution of fragment sizes as 



time elapses. In a recent study [23|, |2J|, using Molec- 
ular Dynamics (MD) simulation with a Langevin ther- 
mostat, we observed a rather complex interplay between 
the polymer chain dynamics and the resulting bond rup- 
ture probability. A major role in this was attributed to 
the one-dimensional (ID) topological connectivity of the 
linear polymer. Significant change in rupture kinetics is 
observed while polymer architecture is tuned as in the 
case of thermolysis of adsorbed bottle-brushes [1^, [2^ . 

Understanding the interplay between elastic and frac- 
ture properties is even more challenging and important 
in the case of 2D polymerized networks (elastic-brittle 
sheets). A prominent example of biological microstruc- 
ture is spectrin, the red blood cell membrane skeleton, 
which reinforces the cytoplasmic face of the membrane. 
In erythrocytes, the membrane skeleton enables it to un- 
dergo large extensional deformations while maintaining 
the structural integrity of the membrane. A number of 
studies, based on continuum- 27 1, percolation- [28l - [30j . 
or molecular level ISll S2l considerations of the mechan- 
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ical breakdown of this network, modeled as a triangu- 
lar lattice of spectrin tetramers, have been reported so 
far. Another example concerns the thermal stability of 
isolated graphene nanoflakes. It has been investigated 
recently by Barnard and Snook [s^ using ab initio quan- 
tum mechanical techniques whereby it was noted that 
the problems "has been overlooked by most computa- 
tional and theoretical studies". Many of these studies 
can be viewed in a broader context as part of the prob- 
lem of thermal decomposition of gels [3J], epoxy resins 
^35*, "36*1 and other 3D networks, studied both experimen- 
tally ;34- 36j|, and by means of simulations [37] as in the 
case of Poly-dimethylsiloxane (PDMS). In most of these 
cases, however, mainly a stability analysis is carried out 
whereas still little is known regarding the collective mech- 
anism of degradation, the dependence of rupture time on 
system size, as well as the decomposition kinetics, es- 
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pecially as far as (2D) polymer network sheets are con- 
cerned. It is also interesting from the standpoint of ba- 
sic physics to compare the degradation process to the 
one taking place in linear polymers where considerable 
progress has been achieved recently (23j . Therefore, in 
the present work we extend our investigations to the case 
of (2D) polymer network sheets, embedded in 3D-space, 
and study as a generic example the thermal degradation 
of a suspended membrane with honeycomb orientation. 

The paper is organized as follows: after a brief intro- 
duction, we sketch our model in Sec. HT] In Sec. IIIII we 
present our simulation results, that is, the distribution 
of bond scission rates over the membrane surface, the 
Mean First Breakage Time (MFBT) of a bond depending 
on membrane size and temperature, the distribution of 
recombination events - Sec. IIII Al and the temporal evo- 
lution of the fragmentation process - Sec. IIII Bl We also 
develop a theoretical scheme based on a set of f*— order 
kinetic differential equations, describing the time varia- 
tion of the number of network nodes, connected by a par- 
ticular number of bonds to neighboring nodes - Sec. IIII CI 
We demonstrate that the analytical solution of such sys- 
tem provides a faithful description of the fragmentation 
kinetics. In Sec. IIVI we conclude with a brief discussion 
of our main findings. 



II. METHODS 
A. Model 

We study a coarse-grained model of honeycomb mem- 
brane embedded in three-dimensional (3D) space. In this 
investigation we consider generally symmetric hexagonal 
membranes {flakes) (Fig. [T|). In a very few cases we 
also discuss fracture of a ribbon shape membranes. The 
membrane flake consists of N spherical particles (beads, 
monomers) of diameter a connected in a honeycomb lat- 
tice structure whereby each monomer is bonded with 
three nearest-neighbors except of the monomers on the 
membrane edges which have only two bonds (see Fig. [T] 
[upper panel]). The total number of monomers N in 
such a membrane is = 6L^ where by L we denote the 
number of monomers (or hexagon cells) on the edge of 
the membrane (i.e., L characterizes the linear size of the 
membrane). There are altogether Nhonds = (3A^ — 6L)/2 
bonds in the membrane. 

We find it appropriate to divide the two-dimensional 
membrane network so that all the beads and bonds are 
distributed into different subgroups presented by concen- 
tric "circles" with consecutive numbers (see Fig. [1] [upper 
panel] ) proportional to the radial distance from the mem- 
brane center. To odd circle numbers thus belong beads 
and bonds that are nearly tangential to the circle. Even 
circles contain no beads and only radially oriented bonds 
(shown to cross the circle in Fig. [TJ The total number of 
circles C in the membrane of linear size L is found to be 
C = {2L — 1). We use this example of dividing the beads 
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FIG. 1: [upper panel] A model of a membrane with honey- 
comb structure that contains a total of A' = 54 beads and has 
linear size L = 3 (L is the number of beads or hexagonal cells 
on the edge of the membrane), [lower panel] An example of 
subdivision of beads and bonds, composing a membrane with 
L — 3, into subgroups ("circles"). The total number of circles 
C in the membrane of linear size L is C = 2L — 1. 



and the bonds composing the membrane in order to rep- 
resent our simulation results in appropriate way which 
relates them to their relative proximity to membrane's 
periphery. 



B. Potentials 

The nearest-neighbors in the membrane are connected 
to each other by " breakable bonds" described by a Morse 
potential, where r is a distance between the monomers, 

Uuir) = eAf {1 - cxp[-Q;(r - rnun)]V (1) 

a = 1 is a constant that determines bond elasticity and 
?'min = 1 is the equilibrium bond length. The dissoci- 
ation energy of a given bond is em = 1, measured in 
units of fcsT, where fc^ denotes the Boltzmann constant 
and T is the temperature. The minimum of this po- 
tential occurs at r = rmin, C^Morse(''inin) = 0. The 

maximal restoring force of the Morse potential, /max = 
^dUyi/dr = aeul'^i is reached at the inflection point. 
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r-infiex = rmin + a ^ ln(2) « 2.69. This force /max de- 
termines the maximal tensile strength of the membrane. 
Stretching of the bond beyond ri„fie^ means potentially 
a scission of that bond as far as the restoring force de- 
clines rapidly with bond length even though there is a 
chance for a recombination. Therefore we take as a crite- 
rion for breaking a bond its expansion to r/i = 5 beyond 
which practically a bond recombinantion is ruled out. 
Since ?7m(0) « 2.95, the Morse potential, Eq. ([1]), is only 
weakly repulsive and beads could partially penetrate one 
another at r < Tmin- Therefore, in order to allow prop- 
erly for the excluded volume interactions between bonded 
monomers, we take the bond potential as a sum of J7m('') 
and the so called Weeks-Chandler- Anderson (WCA) po- 
tential, C/wca(?'), (i.e., the shifted and truncated repul- 
sive branch of the Lennard- Jones potential). 




for r < 2i/V 
for r > 2i/V 



(2) 



with parameter e = 1 and monomer diameter a = 
2-1/6 « 0.89 so that the minimum of the WCA po- 
tential coincides with the minimum of the Morse po- 
tential. Thus, the length scale is set by the parameter 
rniin = 2^/6(7 = 1. The nonbonded interactions between 
monomers are taken into account by means of the WCA 
potential, Eq. The nonbonded interactions in our 

model correspond to good solvent conditions whereas the 
bonded interactions make the bonds in our model break- 
able so they undergo scission at sufficiently high T. 

We have been anxious to emphasize the common fea- 
tures of failure in materials with similar architecture but 
largely varying elasticity properties, e.g., from 1000 GPa 
graphene's Young modulus Q compared to 4 x lO^-^ GPa 
for spectrin [s^l- Putting the value of a Kuhn seg- 
ment to cr = 1.44 A and taking the thermal energy 
knT = 4 X 10-21 J at T = 300 K, we get from our sim- 



ulation [4l[ a Young modulus 0.03 GPa which is ranged 
between typical values for rubber-like materials 0.01- 
0.1 GPa. 



is the three dimensional vector of random force acting 
on particle i. The random force '^i, which represents 
the incessant collision of the monomers with the solvent 
molecules, satisfies the fluctuation-dissipation theorem 
{Ria{t)Rjj3{t')) = 2"fkBT5ij5ai35{t — t') where the symbol 
(. . .) denotes an equilibrium average and the greek-letter 
subscripts refer to the x, y ot z components. The friction 
coefficient 7 of the Langevin thermostat is generally set to 
7 = 0.25 and only in very few cases to 10. The integration 
step is 0.002 time units (t.u.) and the time is measured in 
units of Tinin \/m./eM- We emphasize at this point that in 
our coarse-grained modeling no explicit solvent particles 
are included. In this work the velocity- Verlet algorithm 
is used to integrate the equations of motion. 

Our MD simulations are carried out in the following 
order. First, we prepare an equilibrated membrane con- 
formation, starting with a fully flat configuration shown 
schematically in Fig. [TJ where each bead in the network 
is separated by a distance rmin = 1 equal to the equi- 
librium separation of the bond potential {Um + t/wCA) 
[see Eq. ([T]), and Eq. ©]. Then we start the simulation 
with this prepared conformation and let the membrane 
equilibrate in the heat bath for a very long period of time 
(« 10*" integration steps) at a temperature low enough 
that the membrane stays intact. Fig. [5J (this equilibra- 
tion is done in order to prepare different starting confor- 
mations for each simulation). Then the temperature is 
raised to the working temperature and the membrane is 
equilibrated for 20 MD t.u. (10^ integration steps) which 
interval was found as sufficient to establish equipartition 
(uniform distribution of the temperature throughout the 
membrane) . The time then is set to zero and we continue 




C. MD algorithm 

In our MD simulation we use Langevin dynamics, 
which describes the Brownian motion of a set of interact- 
ing particles whereby the action of the solvent is split into 
slowly evolving viscous (frictional) force and a rapidly 
fluctuating stochastic (random) force. The Langevin 
equation of motion is the following: 



(3) 



where m denotes the mass of the particles which is 
set to TO = 1, iti is the velocity of particle i, ~^ i = 
C^M + ^^wcA)i is the conservative force which is a 



sum of all forces exerted on particle i by other parti- 
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cles in the system, 7 is the friction coefficient and Ri 



FIG. 2: A snapshot of a typical conformation of an intact 
membrane with L = 30 containing 5400 monomers after equi- 
libration. Characteristic ripples are seen to cross the surface. 



the simulation of this membrane conformation to exam- 
ine the thermal scission of the bonds. We measure the 
elapsed time r until the flrst bond rupture occurs and 
repeat the above procedure for a large number of events 
(lO'^-lO'^) so as to sample the stochastic nature of rup- 
ture and to determine the mean (r) which we refer to as 
the Mean First Breakage Time. In the course of simula- 
tion we also calculate properties such as the probability 



distribution of breaking bonds regarding their position 
in the membrane (a rupture probabihty histogram), the 
probabiUty distribution of the first breakage time W{t), 
the mean extension of the bonds with respect to the con- 
secutive circle number in the membrane, as well as other 
quantities of interest. 

Since in the problem of thermal degradation there is 
no external force acting on the membrane edges, a well- 
defined activation barrier for a bond scission is actually 
missing, in contrast to the case of applied tensile force. 
Therefore, a definition of an unambiguous criterion for 
bond breakage is not self-evident. Moreover, depend- 
ing on the degree of stretching, bonds may break and 
then recombine again. Therefore, in our numeric exper- 
iments we use a sufficiently large value for critical ex- 
tension of the bonds, = 5r,iii,i, which is defined as a 
threshold to a broken state of the bond. This convention 
is based on our checks that the probability for recom- 
bination (self-healing) of bonds, stretched beyond rh, is 
sufficiently small, < 10^^ . 

Also we examine the course of the degradation kinet- 
ics: at periodic intervals in separate simulations runs we 
analyze the size distribution of fragments (clusters) of the 
initial membrane and establish the time-dependent prob- 
ability distribution function of fragment sizes, P(n,i), as 
time elapses after the onset of the thermal degradation 
process. This also yields the time evolution of the mean 
fragment size, N{t) = n{t)P{n,t)dn. We perform 
the statistical averaging of fragment sizes by developing 
an appropriate for the system fast cluster counting algo- 
rithm. 

The Molecular Dynamics calculations were carried out 
using a cluster counting algorithm based on an ad hoc 
implementation of the Hoshen-Kopelman program [ssj ]. 
fn a subsequent paper of Al-Futaisi and Patzek [s^, the 
Hoshen-Kopelman algorithm for cluster labeling has been 
extended to non-lattice environments where network el- 
ements (sites or bonds) are placed at random points in 
space. Following Al-Futaisi and Patzek [s^, we devel- 
oped our simplified version of this algorithm which con- 
cerns only clusters of nodes in a network with honeycomb 
structure (which is our model shown in Fig. [5]). In our im- 
plementation we assume that all nodes in the network are 
occupied, but some links can be broken when we study 
thermal degradation process of the network. The net- 
work information is stored in two arrays. The first one is 
related to the connectivity of the nodes and the second 
to the state of links (intact or broken link) 



III. RESULTS AND DISCUSSION 

A. Bond scission 

In Fig. [3^ we show the distribution of bond scission 
rates among all bonds of the honeycomb membrane for 
flakes (with zig-zag pattern on the periphery) of several 
sizes L = 10, 15, 20. Somewhat surprisingly, one finds 




Consecutive circle number 

FIG. 3: (a) Rupture probability histograms for thermally in- 
duced scission events in flexible honeycomb flake of differ- 
ent linear size L as shown in the legend, (b) Rupture his- 
togram for a ribbon-like square honeycomb membrane with 
496 nodes, (c) Rupture histogram for a membrane flake with 
beads tethered at the rim, L = 10. Parameters of the heat 
bath are T = 0.1 and 7 = 0.25. (d) Probability distribution 
of breakage events as a function of consecutive circle num- 
ber for a membrane flake with A'^ = 600 and two different 
friction coeflicients 7 = 0.25 and 10.0. Here T = 0.1. The 
inset shows estimated values of Lyapunov exponents A vs T 
for beads located in the rim/bulk of membrane as indicate. 
Here iV = 5400, 7 = 0.25. 



the overwhelming fraction of bond breaking occurs at the 



5 



outer-most rim of the membrane where monomers are 
bound by only two bonds to the rest of the sheet. We 
have also sampled the rupture histograms for ribbon-like 
square membranes, Fig. [5]d. Interestingly, we observe no 
difference between the rupture rates of zig-zag and arm- 
chair edges whereas the bonds in the four corners of such 
membrane expectedly break more frequently. The differ- 
ence in the relative stability of the bonds becomes clearly 
evident in Fig.[3ji where the frequency of periphery bonds 
appears nearly two orders of magnitude larger when com- 
pared to bonds in the 'bulk' of the membrane where each 
monomer (node) is connected by three bonds to its neigh- 
bors. One may therefore conclude that a moderate in- 
crease in the coordination number of the nodes (by only 
33% regarding the maximum coordination of a node) 
leads to a major stabilization of the supporting bonds 
and much stronger resistance to fracture. Our additional 
simulation in the strongly damped regime for 7 = 10 
indicates no qualitative changes compared to 7 — 0.25 
except an absolute overall increase of the rupture times 
which is natural for a more viscous environment. 

Note that the question of where and which bonds pre- 
dominantly break is by no means trivial. For example, 
in the case of linear polymer chain thermal decompo- 
sition the rate of bond rupture is least at both chain 
ends although the end monomers, in contrast to those 
inside the chain, are bound by a single bond only [23j . 
We demonstrate in Fig. ^ that this interesting feature 
holds also for the honeycomb membrane flake, provided 
the rim is clamped and left immobile during the simula- 
tion. Evidently, the highest frequency of bond scissions 
grows towards the membrane center. Similar to the case 
of polymer chains with fixed ends [2^, . 

In order to provide deeper insight into the mechanism 
of temperature-induced bond breaking, in the inset to 
Fig. [3JI we present the temperature variation of Lya- 
punov's exponent A for membrane nodes located in the 
bulk and in the rim of the sheet. Evidently, beyond a 
cross-over temperature T « 0.05 one observes a signifi- 
cant growth of Ari„ as compared to Abuik- This indicates 
that the trajectories of nodes at the membrane periphery 
attain much faster chaotic features at higher temperature 
than those of the bulk nodes. Moreover, we should note 
that beads in the vortices have values of A which exceed 
those in the rim by about 5%. Therefore this analysis 
of trajectory stability at characteristic locations in the 
membrane clearly demonstrate that bond rupture is in- 
duced by intermittent motion of the respective nodes. 

The variation of the MFBT r of a bond with mem- 
brane size N during thermolysis for both hexagonal and 
square shapes of the 2D sheet is displayed in Fig.|4^. Ev- 
idently, one observes for r a well pronounced power law 
behavior, r c>c 7V~'^ with an exponent /3 « 0.50 ± 0.03. 
It turns out that the scaling exponent /? remains insensi- 
tive to changes in the geometric shape of the membrane 
sheet. This value of /3 might appear somewhat surpris- 
ingly to deviate from the expected exponent of unity, 
given that in the absence of external force all bonds are 



supposed to break completely at random so that the total 
probability for a bond scission (i.e., the chance that any 
bond might break within a time interval) is additive and 
should be, therefore, proportional to the total number of 
available bonds, N^onds = (3A^ — 6i)/2. As suggested 
by Fig. 131 however, predominantly only periphery bonds 
are found to undergo scission during thermal degrada- 
tion. The number of periphery bonds goes roughly as 
cx which agrees well with the observed value /3 w 0.5 
and provides a plausible interpretation of the simulation 
result. Fig. From the inset in Fig. one may verify 
that the bond scission displays an Arrhenian dependence 
on inverse temperature, r cx exp^AEb / ksT) , with a slope 
AEb ~ 1. This slope suggests a dissociation energy AEb 
of the order of the potential well depth of the Morse in- 
teraction, Eq. ([1} where sm = 1-0. In our model we deal 
typically with Eb/iksT) w 10 which at 300 K and typi- 
cal bond length r„i„ w 0.144 nm, corresponds to ultimate 
tensile stress ~ 0.6 GPa. This is a reasonable value for 
our membrane which is considerably softer than graphene 
with - 100 GPa and is ranged between typical values 
for rubber materials 0.01-0.1 GPa. 

The probability distribution function (PDF) of MFBT 
W{t), i.e., the PDF of the time interval before the first 
breakage event in the membrane takes place, is shown 
in Fig. |3Jd for T = 0.1. Evidently, one has W{t) cx 
exp(— i/r) with a sharp maximum close to < « 0. At 
even lower temperature one might expect this maximum 
to become more pronounced, suggesting thus a Poisson 
distribution for the W{t). As far as r tends to grow ex- 
ponentially fast with decreasing T - see inset in Fig. 
collecting statistics in this temperature range becomes 
difficult. 

We also analyze the fracture of bonds with respect to 
their possible recombination. The procedure of sampling 
bond recombination events is the following. Once one of 
the bonds in the membrane is stretched to the position of 
the Morse potential infiexion point r^^u^^ w 2.69 (corre- 
sponding to maximal tensile strength of a given bond), we 
start to monitor its further expansion for 10^ integration 
steps and record its maximal expansion h — r — r^^a^^. 
If the bond subsequently shrinks to r < r^^a^^, this 
is counted as recombination event. Simultaneously we 
record the time t needed to go back below the expansion 
''intiox- We then compute distributions Qh{h) of bond ex- 
pansions beyond r^^a^^ and recombination times Ph(t). 
Thus Qh{h) yields the probability of bond overstretching 
to distance h ~ r — r^^^fi^^ during the recombination event. 
As indicated by the distribution Ph {t) shown in Fig. [SJ 
the interval 20 t.u. exceeds more than five times the max- 
imal life of recombination which guarantees that all such 
events are properly taken into account. Generally we find 
that a recombination of bonds takes place rather seldom 
(roughly once per simulation run of average length 130 
t.u., as given by the MFBT). 

It has been mentioned in Sec. |TT] that in the absence 
of external force acting on the membrane the criterion 
for a bond to be considered broken is not unambigu- 
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FIG. 4: (a) Mean first breakage time (r) vs. number of beads 
A'^ for two different membrane siiapes: a fiexagonal flake and 
a square ribbon. Solid line represents a fit by power law with 
an exponent /3 — 0.5 in both cases. The inset shows the varia- 
tion of (r) with inverse temperature 1/T for a flake membrane 
with A*' — 294 particles. The fitting line yields an Arrhenian 
relationship, (r) oc exp{AEi, / kBT) with dissociation energy 
barrier AEb ~ 0.95. (b) Breakage time probability distribu- 
tions W{t) for thermally induced breakage of a flake mem- 
brane with A' = 600 particles at T = 0.1. Symbols denote 
result of the simulation and full line stands for a fitting func- 
tion W{t) oc exp{-t/T) with r = 128.2. 



ous. Adopting a rather large critical bond extension of 
rh — 5r„i„ practically rules out the probability of of sub- 
sequent recombination of such bonds, as indicated by 
Qh{h) oc exp(— 3.2/i) in [5] (right inset). Indeed, if one 
adopts ri„fie^ as a threshold for rupture than Qh{h) sug- 
gests that bond expansions larger than r^^a^^ -f 1.1 ~ 3.79 
practically never happen. The PDF P{t) of times elapsed 
before self-healing is also shown in Fig.[S]where it appears 
as a Poisson distribution, P{t) « texp(— 0.3^). Not sur- 
prisingly, most of the recombination occur in the periph- 
ery of the membrane - left inset in Fig. [SJ 

In Fig. [S] we show the distribution of broken bonds 
at different times after the onset of thermal degrada- 
tion. We would like to note that the regular hexagonal 
flakes shown in Fig. |6] serve only to indicate schemati- 
cally the positions of both broken and still intact bonds 
and by no means represent the actual conformation of 
the membrane. Two different temperatures, T = 0.10 




FIG. 5: Probability distribution Ph{t) of times (impulses), 
and Qh{h) of bond overstretching lengths h (circles, right 
panel of inset) before a recombination event in a membrane 
with N = 600, T = 0.1, 7 = 0.25 takes places. The re- 
combination times probability distribution Ph{t) is fitted by 
Poisson distribution (blue line). The probability for a bond 
stretching a distance h beyond ri„ticx is described by Qh{h) in- 
dicating an exponential decay - red line. The left inset shows 
healing probability Rh vs. consecutive circle number, demon- 
strating which part of the membrane undergoes healing most 
frequently. 



and T = 0.15 are studied. One can readily verify from 
Fig.[6^-d that a.tT — 0.10 the degradation process starts 
from the rim of the network sheet and then proceeds in- 
wards, as noted earlier by Meakin (4^. In contrast, at 
50% higher temperature, i.e., at T = 0.15, bonds break 
everywhere in the network sheet - cf. Fig. [6|d, f and 
Fig. [nt, g. Such difference in the bond scission mech- 
anism at different temperatures has been observed [i^ 
before. It appears that at T = 0.1 the membrane periph- 
ery undergoes stronger oscillations than the membrane 
bulk which lead to bond scission at the rim while in the 
inner part of the sheet the monomers mutually block 
each other and the bonds remain largely intact. This 
explanation is supported by the measured values of Lya- 
punov exponents (see the inset of Fig. [3] and discussion 
in SecIiniA). On the other hand a 'hot' i.e. at T = 0.15 
membrane undergoes much stronger agitation as a whole 
and as a result bonds break all over the network sheet. 
Moreover, at lower T when the thermal agitation of the 
membrane is weaker, one may correlate the probability 
of bond scission with the average strain in the bonds - 
Fig. [71 Evidently, the average squared bond length (Z^) 
increases steadily and becomes nearly 8% larger for the 
bonds sitting on the last ring of the membrane whereas no 
difference in the mean kinetic energy between peripheral 
and bulk nodes is detected - Fig. [71 Note that the kinetic 
energy (that is, the temperature T) is thereby uniformly 
distributed along the sheet - see inset in Fig. [71 
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FIG. 6: Thermal breakage of bonds in a membrane made of 
— 600 particles at different time moments: a) t — 10^, b) 
7 • 10^ c) 5 • 10^ d) 3 • lO'^ and e) 50, f) 100, g) 250, h) 500. 
Broken bonds are marked by blue (T = 0.1) or red (T = 0.15) 
color, depending on the temperature of a heat bath, while gray 
color corresponds to intact bonds. 



B. Temporal evolution of the fragmentation 
process 

After the onset of the thermal decomposition process 
the membrane flake disintegrates into smaller clusters 
(fragments) of size n whose mean size (or average molec- 
ular weight) N{t) decreases steadily with time. N{t) 
is easily accessible experimentally, therefore we give in 
Fig. the course of its temporal evolution, observed in 
our computer experiment. Using an ad hoc cluster count- 
ing program in the course of the MD-simulation we sam- 
ple the probability distribution of fragment sizes, P{n, t), 
so that the first moment N{t) = J n{t)P{n, t)dn gives the 
cluster mean size N(t). Thus, for a given time moment 




Consecutive circle number 



FIG. 7: Variation of mean squared bond length {l^) vs. dis- 
tance from membrane center (consecutive circle number) for 
thermolysis of membrane composed of A*' = 5400 beads. The 
inset displays mean kinetic energy (Ek) (per monomer) as 
a function of circle number. Red dashed line represents the 
level of energy which corresponds to equipartition theorem. 
Parameters of thermostat are T = 0.1 and 7 = 0.25. 



t we average data over more than 10'^ independent runs, 
each starting from a different initial conformation of the 
honeycomb membrane. In Fig.[5)D we show the time vari- 
ation of the ensuing PDF P{n, t) whereby the system is 
seen to start with a single sharp peak at t = when 
the membrane is still intact. As time goes by, P{n, t) 
becomes bimodal, the maximum of the distribution is 
seen to shift to smaller values of cluster size whereas an 
accumulation of fragments of size 1 or 2 is observed to 
contribute to a second peak at n « 1. Eventually, as 
t — 7> 00, the PDF P{n, t) settles to a shape with a single 
sharp peak (a (5-function) at n « 1 (not shown here). 

One can readily see from Fig. |51i that the quantity 
1 — N~^{t) does not immediately follow a straight line 
of decay when plotted in semi-logarithmic coordinates, 
rather, such a decay is observed after an initial period 
of slower decline. This effect is due to averaging over 
many realizations of the fragmentation process. In each 
run the degradation of bonds starts earlier or later at a 
particular time r (the Mean First Breakage Time) that is 
distributed according to W{t) - cf. Fig.|3}D. As a result a 
clear cut exponential course of 1—N~^{t) is only observed 
in the late stages of fragmentation. Such behavior is 
found independently of the membrane size - Fig. [5^. 

In addition, one could expect that the fragmentation 
process is not governed by a single rate constant in 
a presumably 1^*— order chemical reaction even though 
the bonds that undergo rupture are chemically identical. 
Therefore, from the temporal mean cluster size behav- 
ior, presented in Fig. one may conclude that even in 
the case of a homogeneous membrane the thermal degra- 
dation process is more adequately described by several 
reaction constants which govern the dissociation of dif- 
ferent groups of bonds. In the next section we suggest 
a simple model of reaction kinetics which takes into ac- 
count this conjecture. 
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several I'^'-order chemical reactions. The nodes of the 
membrane can be subdivided into several groups, de- 
pending on the number of intact bonds that connect them 
to neighboring nodes. In the case of a honeycomb mem- 
brane one may distinguish four such groups, and denote 
the instantaneous number of such nodes (monomers) by 
no, ni, 7i2 and whenever 0, 1, 2, or 3 intact bonds 
exist around such a node. If self-recombination of bonds 
is ignored, which is reasonable in view of large value of 
the threshold, and simultaneous scission of two and more 
bonds is disregarded as hardly probable, one can write 
down a strongly simplified set of l**-order kinetic dif- 
ferential equations (DE) that describes the evolution of 
rto(i), rii(t), n2{t) and 713(1) with time: 
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ni{t) = ~kini{t) + k2n2{t) 
n2{t) = -k2n2{t) + ksnsit) 



(4) 



Thus, for example, the number of double-bonded nodes 
decreases as one of the two bonds breaks, however, if a 
bond breaks around a node with triple coordination this 
would increase the population of double-bonded nodes. 
Note, that the total number of nodes of all kinds remains 
thereby conserved. 



no{t) + ni{t) + n2{t) + n^it) = N ^ GL"^ 



(5) 



FIG. 8: (a) Semi-log plot of time variation of mean fragment 
size N{t) for membranes made of different size as indicated. 
Symbols represent simulation results whereas a red line stands 
for the fitting function 1 — 1/N(t) oc exp{ — kt). The kinetic 
constant k = 4.3~^s^^. In the inset the same is shown in 
normal coordinates, (b) Probability distribution of fragment 
sizes P{n, t) at different times t (in MD time units) after be- 
ginning of the thermal degradation process for a membrane 
with A'^ = 294. Parameters of a heat bath are T = 0.12 and 
7 = 0.25. 



C. Reaction kinetics 

One may try to reproduce the bond scission kinetics 
during in the course of thermal fragmentation by a set of 



If the degradation process starts with an intact mem- 
brane conformation at t = (no broken bonds exist), 
one can fix the initial conditions as no(0) = 0, ni(0) = 
0,712(0) = 6L and 713(0) =N-n2{0) = 6L(L-1). Thus, 
initially one has only ^2(0) = 6L double-bonded nodes 
at the membrane periphery along with 713(6) — 6L{L—1) 
triple-bonded nodes in the bulk of the flake. Then, as the 
decomposition process develops, nodes from a given class 
r — 1, 2, 3 will transform into a lower class r = 0, 1, 2 
ones (we neglect hereby the simultaneous scission of more 
than one bond of a node as highly improbable event). 
Eventually, at i 00 the fragmentation process ends 
and one expects 7ii(oo) = 712(00) = 713(00) = and 
7Zo(oo) = N. 

One may solve analytically the system of l^*-order DE 
Eqs. (11]) to a set of functions ni{t) with < z < 3: 



no{t) = 6L^ - ni{t) - n2{t) ~ ns{t) 

n,{i) = , rr [^l(fe-fc3)(e-'=^*-e-'=^*)+fe^3(e-'■^*-e- 
(fcl - fc2)(fci - fc3)(fc2 - fc3) 
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(6) 



r 



where the rate constants fci, ^2 and fcs are stih to be 
determined, for example, by comparison with simulation 
data. As one may readily verify from Fig. [3J our simu- 
lation data suggest that a bond to a triple-bonded node 
at T = 0.12 is much more stable (by about two orders 
of magnitude) than a bond at the flake periphery where 
each node is connected by only two bonds to the net- 
work. Thus the system of Eqs. (6) can be tested for a set 
of reaction constants ki ^ k2 > k^ directly by means of 
our Molecular Dynamics computer experiment. 




FIG. 9: (a) Variation of the number of nodes nm{t) with 
m = 0-3 intact bonds with elapsed dimensionless time t/riT) 
for a membrane with A'^ = 54 and T = 0.10. Here r{T) 
denotes the characteristic time of degradation at the respec- 
tive temperature T. The ni{t) (monomers bound by a single 
bond) are shown by shaded area, (b) The same as in (a) for 
the same membrane size A'^ = 54 according to the analytic 
result Eqs. (6). The values of the kinetics rate constants are 
ki = 20.0, fca = 0.007, fcg = 0.005. 



Indeed, one finds by comparing the simulation result. 
Fig. for a membrane of size A^ = 54atr = 0.10, and 
the analytical solution, Eqs. (6), Fig. IHb, that the ob- 



served kinetics agree qualitatively, provided one allows 
for the absence of fluctuation in Eqs. (4) (i.e., for the 
disarray and averaging of the MFBT r in the simulation 
data). We should like to point out here that the values 
of the rate constants ki, k2, k^ are not best fit values. 
Because of the additional effect of averaging discussed in 
Sec. nil Bl and shown in Fig. [SI a straight fitting procedure 
with three parameters ki,k2, ks would be both inefficient 
and hardly successful. Therefore we tried different com- 
binations of values for fci, fc2, ks subject to the condition 
that ki k2 > kz- Even though the general qualita- 
tive shapes turn out to be rather sensitive to the values 
of ki,k2, k^ (being capable of reproducing, e.g., the exis- 
tence of a common intersection point of no(^),f^2(^) and 
n^{t) - see Fig. [5^ - for particular choice of parameters), 
one finds thus easily a combination which qualitatively 
matches well the simulation data shown in Fig. 

One may conclude therefore that the simplified set of 
1***— order DE Eqs. (HI) captures qualitatively the main 
features of the fragmentation kinetics and the principal 
mechanism at work is a combination of few 1"*— order 
chemical reactions of bond scission. Nonetheless, it is 
conceivable to expect that for a full quantitative descrip- 
tion of the thermal degradation process the set of kinetic 
equations, Eqs. (4), should be extended by few additional 
reactions: ni{t) — n2{t), n2{t) — n^{t) that may in 
principle also take place (with respective rate constants). 
One can then still derive an analytical solution of the ex- 
tended set of DE that describes the full kinetics of frag- 
mentation and try to fit the ensuing rate constants to the 
simulation data. In view of the growing number of fit pa- 
rameters, however, a detailed analysis of such system is 
beyond the scope of the present work and should be left 
for future work. 



IV. CONCLUSION 

In the present investigation we use Langevin Molecu- 
lar Dynamics simulation and also solve a set of 1***— order 
kinetic DE so as to model the process of thermal destruc- 
tion of a polymerized membrane sheet with honeycomb 
structure. Results of our work can be treated as generally 
applicable due to the fact that a two-dimensional regu- 
lar lattice can only be created in few ways: honeycomb, 
triangular and square lattice with second neighbor inter- 
actions because of the zero-shear modulus (apart from 
exotic cases like quasi-crystalline, kagome, etc. lattices 
of little relevance) . The differences in lattice coordination 
among the first three periodic networks induce only quan- 
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titativc rcnormalization of the Young modulus. However 
they don't change the overall elastic behavior. Our find- 
ings regarding the most salient features of thermolysis 
in an elastic brittle honeycomb network sheet subject to 
sufficiently high temperature can be summarized as fol- 
lows: 

• The probability of bond scission is highest at the 
periphery of the membrane sheet where nodes are 
connected by two bonds only. At higher tempera- 
ture, however, the whole sheet undergoes fragmen- 
tation whereby also bonds in the bulk rupture. 

• The mean time r until a bond undergoes scission 
event declines with the number of nodes N (with 
membrane size) by a power law as r oc A^^*^-^ in- 
dependently of the geometric shape of membrane 
sheet. The times of bond scission are exponentially 
distributed, W{t) oc exp(— ^/(t)). 

• Bond recombination (self-healing) occurs seldom 
during thermal degradation, and the measured 
recombination times follow Poisson distribution 
whereas an extra-stretching of bonds (beyond the 
point of maximal tensile strength) before a recom- 
bination takes place is exponentially improbable. 

• the fragmentation kinetics is determined by l***- 

order reactions between network nodes with dif- 
ferent number of intact bonds and follows simple 



exponential decay at late times. 

• A set of l**-order kinetic differential equations, de- 
scribing the process of fragmentation, can be es- 
tablished and solved analytically. One finds re- 
sults in qualitative agreement with those from com- 
puter experiment providing thereby a deeper in- 
sight into the mechanism of thermal degradation of 
two-dimensional honeycomb networks. 

In view of the presented results, it should be clear that 
more work is needed (e.g., regarding the effects of ran- 
domness in brittle networks) until a full understanding of 
the process of thermal degradation in polymerized 2D- 
membranes is reached. 
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